#include "CompressedColumnStorage.h"
#include <stdio.h>

CompressedColumnStorage::CompressedColumnStorage(const DMatrix& _M){
    
    m_dim = _M.getCols();
    int non_zeros = m_dim*(m_dim+1)/2;
    
    
    //total of non_zeros elements:
    for(int j=0; j<m_dim; j++){
        for(int i=0; i<j; i++){
            if(_M(i,j) == 0){
                non_zeros--;
            }
        }
    }
    
    
    row_index.resize(non_zeros);
    col_ptr.resize(m_dim + 1);
    value.resize(non_zeros);
    
    //storage method:
    int index=0;
    col_ptr(0) = 0;
    for(int j=0; j<m_dim; j++){
        for(int i=0; i<=j; i++){
            if(_M(i,j) != 0){
                value(index) = _M(i,j);
                row_index(index++) = i;
            }
            
        }
        col_ptr(j+1) = index; 
    }
    
}

CompressedColumnStorage::~CompressedColumnStorage()
{
 
 
}

double CompressedColumnStorage::getElement(int i, int j){

    // Find where the row starts
    int columnStart = col_ptr(j);

    // Find where the next row starts
    int columnEnd = col_ptr(j+1);

    // Scan the column indices of entries in row i
    for(int k = columnStart; k < columnEnd; k++)
    {
            // if j is the column index of a non-zero, then return it
            if(row_index(k) == i)
                    return value(k);

            // if we have passed j, then entry (i, j) must be a zero
            if(row_index(k) > i)
                    return 0;
    }

    // if we have reached the end of the non-zeroes without
    // find j, it must be the index of a trailing 0
    return 0;
}

//print the matrix in CompressedRowStorage format
void CompressedColumnStorage::print(){
    cout << "Values vector:\n";
    for(int i=0; i<value.getLength(); i++){
             cout << value(i) << "\n";
    }
    cout << "\nRow Index Vector:\n";
    for(int i=0; i<row_index.getLength(); i++){
            cout <<row_index(i) << "\n";
    }
    cout << "\nColumn Pointer Vector:\n";
    for(int i=0; i<col_ptr.getLength(); i++){
            cout << col_ptr(i) << "\n";
    }
}

DMatrix CompressedColumnStorage::pre_multiply(const DMatrix &_M){
    DMatrix result_matrix(m_dim,_M.getCols());
    result_matrix.zero();
    for(int j=0; j<m_dim;j++ ){
        for(int k=0; k<_M.getCols(); k++){
            for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
                result_matrix(j,k) += value(i)*_M(k,row_index(i));
                result_matrix(row_index(i),k) += value(i)*_M(k,j);
            }
            result_matrix(j,k) += value(col_ptr(j+1)-1)*_M(k,j);
        }
    }
    return result_matrix;
}

DMatrix CompressedColumnStorage::post_multiply(const DMatrix& _M) const {
	DMatrix result_matrix(_M.getRows(),m_dim);
	result_matrix.zero();
	for(int j=0; j<m_dim;j++ ){
		for(int k=0; k<_M.getRows(); k++){
			for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
				result_matrix(k,row_index(i)) += value(i)*_M(k,j);
				result_matrix(k,j) += value(i)*_M(k,row_index(i));
			}
			result_matrix(k,j) += value(col_ptr(j+1)-1)*_M(k,j);
		}
	}
	return result_matrix;
}

DVector CompressedColumnStorage::pre_multiply(const DVector& _V){
	DVector result_vector(m_dim);
	result_vector.zero();
	for(int j=0; j<m_dim;j++ ){
			for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
				result_vector(j) += value(i)*_V(row_index(i));
				result_vector(row_index(i)) += value(i)*_V(j);
			}
			result_vector(j) += value(col_ptr(j+1)-1)*_V(j);
	}
	return result_vector;
}


DVector CompressedColumnStorage::post_multiply(const DVector& _V) const{
	DVector result_vector(m_dim);
	result_vector.zero();
	for(int j=0; j<m_dim;j++ ){
		
			for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
				result_vector(row_index(i)) += value(i)*_V(j);
				result_vector(j) += value(i)*_V(row_index(i));
			}
			result_vector(j) += value(col_ptr(j+1)-1)*_V(j);
	}
	return result_vector;
}

DMatrix CompressedColumnStorage::operator*(const DMatrix& _M){
	return this->pre_multiply(_M);
}

DVector CompressedColumnStorage::operator*(const DVector& _V){
	return this->pre_multiply(_V);
}

const DMatrix operator*(const DMatrix& _M, const CompressedColumnStorage& _CS) {
	return _CS.post_multiply(_M);
}


const DVector operator*(const DVector& _V, const CompressedColumnStorage& _CS) {
	return _CS.post_multiply(_V);
}

DMatrix CompressedColumnStorage::toMatrix(){
	DMatrix matrix(m_dim,m_dim);
	matrix.zero();
	for(int j=0; j<m_dim;j++ ){
		for(int i=col_ptr(j); i<col_ptr(j+1)-1; i++){
			matrix(row_index(i),j) = value(i);
			matrix(j,row_index(i)) = value(i);
		}
		matrix(j,j) = value(col_ptr(j+1)-1);
	}
	return matrix;
}
        
        


